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Abstract. 

We propose a general algorithm to enumerate all solutions of a zero-dimensional 
polynomial system with respect to a given cost function. The algorithm is de- 
veloped and is used to study a polynomial system obtained by discretizing the 
steady cavity flow problem in two dimensions. The key technique on which 
our algorithm is based is to solve polynomial optimization problems via sparse 
semidefinite programming relaxations (SDPR) _20J, which has been adopted suc- 
cessfully to solve reaction-diffusion boundary value problems in [13]. The cost 
function to be minimized is derived from discretizing the fluid's kinetic energy. 
The enumeration algorithm's solutions are shown to converge to the minimal 
kinetic energy solutions for SDPR of increasing order. We demonstrate the 
algorithm with SDPR of flrst and second order on polynomial systems for dif- 
ferent scenarios of the cavity flow problem and succeed in deriving the k smallest 
kinetic energy solutions. The question whether these solutions converge to so- 
lutions of the steady cavity flow problem is discussed, and we pose a conjecture 
for the minimal energy solution for increasing Reynolds number. 
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1 Introduction 



The steady cavity flow problem is a simple model of a flow with closed stream- 
lines and is used for examining and validating numerical solution techniques 
in fluid dynamics. Although it has been discussed in several literature of nu- 
merical analysis of fluid mechanics (see, e.g., [S], [3], [7], [1], it is still an 
interesting problem to a number of researchers for a range of Reynolds numbers. 
We are interested in a polynomial system derived from discretizing the steady 
cavity flow problem. This polynomial system, called the discrete cavity flow 
problem, is obtained by discretizing the cavity region, approximating the partial 
differential equation of the two-dimensional cavity flow problem by finite differ- 
ence method, and depends on two parameters, the Reynolds number R and the 
boundary velocity v. 

Our main contribution is an algorithm to enumerate the discrete cavity flow 
problem's solutions with respect to an objective function, that is derived from 
discretizing the cavity flow's kinetic energy function. The key element of the enu- 
meration algorithm is the sparse semidefinite programming relaxation method 
(SDPR) [2D| for solving polynomial optimization problems, whose solution is 
taken as the starting point for Newton's method or sequential quadratic pro- 
gramming. Recently, the SDPR has been successfully adopted to derive numer- 
ical solutions to a class of reaction diffusion equations [13]. In this paper, the 
polynomial optimization problem is the minimization of the discretized kinetic 
energy subjected to the discrete cavity flow problem. We prove that the first 
k solutions provided by the enumeration algorithm converge to the k smallest 
energy solutions of the discrete cavity flow problem, in case that we apply SDPR 
of increasing relaxation order. Furthermore, we demonstrate this algorithm for 
different parameter settings of R and v, and show in some examples that it 
is sufficient to apply SDPR with first or second order, to enumerate accurate 
approximations to the smallest energy solutions. At second, we discuss the 
minimal energy solution's behavior of the discrete steady cavity flow problem in 
case that a finer grid is chosen to discretizc the cavity fiow problem. For small 
Reynolds numbers R standard grid-refining techniques can be applied to extend 
solutions of the polynomial system to finer grids. However the polynomial sys- 
tems for large R and v behave differently and convergence is far more difficult 
to obtain. We examine the polynomial systems for a fixed discretization and 
increasing Reynolds number R. Based on our observations, we conjecture the 
minimal kinetic energy is converging to zero if R tends to infinity. Also, we test 
the performance of SDPR for an alternative finite difference discretization by 
Arakawa [1] of the steady cavity fiow problem. 

The exact formulation of the steady cavity fiow problem in 2 dimensions is 
introduced in section 2. We discuss its boundary conditions and derive the dis- 
crete steady cavity flow problem. In section 3 we study the discrete cavity flow 
problem by Grobner basis method for coarse grid discretizations, in order to 
be able to verify the results derived by our enumeration algorithm later on. In 
section 4 we show how to solve a polynomial optimization problem derived from 
the discrete cavity flow problem via the SDPR method. We present our main 
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contribution, the algorithm to enumerate the discrete cavity flow problems solu- 
tions w.r.t. their kinetic energy, and demonstrate its power for some settings of 
the 2 parameters R and v. In section 5 we discuss the question of how to extend 
the discrete cavity flow problem's minimal energy solution to finer grids and its 
convergence to analytic solutions of the steady cavity flow problem. Finally, we 
examine the discrete steady cavity flow problem for increasing Reynolds number 
R. 



2 A Polynomial System for the Steady Cavity 
Flow Problem 

Let us review the well-known stream function method to solve the Navier-Stokes 
equation (see, e.g., [S], [TB]). The stream function method is a standard method 
to solve the 2 dimensional steady cavity flow problem numerically. 

Let {u{x,y,t),v{x,y,t)) be the velocity of the two dimensional cavity flow 
of an incompressible fluid. It follows from the continuity equation of the in- 
compressible fluid (preservation of the mass) + f| = that there exists a 
fimction ip{x,y,t) such that 



dip 



u. 



(1) 



dip 

dx ' dy 

Put V = (u, u,0). cD ~ rotv is called the vorticity. Since the last coordinate 
of V is 0, cJ can be written as (0, 0, ct;(a;, y, t)). The continuity equation and 
the Navier-Stokes equation (preservation of the momentum) can be written as 
follows in terms of Tp and lo. 

Aijj = ~uj (2) 



duj dip duj dip duj 1 ^ 
dt dy dx dx dy R 



(3) 



Here, A is the Laplace operator and R is the Reynolds number. Let us consider 
the cavity region ABCD with the coordinate A = (0,0), B = (0,-1), C = 
(1,-1), = (1,0). 



A — 



D 



Jc 



The steady cavity flow problem is ([2]) and ([3]) with the steady condition 



^ = and the boundary condition 



u(0, y) = w(x, 0) = u(l, y)^0 on AB, BC, CD 
v{0, y) = v{x, 0) = v{l, y)^0 on AB, BC, CD 
u{x, 1) = s, v{x, 1) = on AD 



(4) 
(5) 
(6) 
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Here s is the velocity of the stream out of the cavity ABCD. 

We devide the square ABCD into the N x N mesh. Put h = 1/N. Let us 
translate these boundary conditions into boundary conditions for ip and ui. It 
follows from (O, ([S]) that the function ijj is constant on the boundaries AB, 
BC, CD, DA. Since ^ is continuous, we suppose that = on. the boundaries. 

The boundary condition for w is a little comphcated; see, e.g., [111 p. 162]. 
We cite this discussion in text books for reader's convinience. Let us consider 
the case of the boundary AD. We take a mesh point M on AD. Let P be the 
mesh point inside the cavity ajacent to M and P' the mirror image of P with 
respect to AD. We supposed that the size of the mesh is h. 



.D 



We denote the value of ■0 at the point P by ■ip{P) or tpp. We have —uj{M) = 
Aip{M) = ipyy ~ ipp-'2-^^M+4'r' _ need to determine the value of ^pp' to get 
an approximate value of lu at M. By using the central difference approximation, 
s = v = |^(M) ~ Then, Vp' ^ 2hs + Tpp. Therefore, we have 

2ijp + 2hs 

j;2 — (7) 

Analogously, we have 

when M is a grid point on AB or BC or CD and P is the adjacent internal 
grid point of M. 

It follows from the discussion above, we obtain the following central differ- 
ence scheme for the steady cavity flow problem. 



li 

+ — {{ipi+ij - Vi-ij){<.^i,j+i - <^i,j-i) (9) 

— {i'ij + l - V'i,J-l){<^i+l,j -<^i-l,j)} 

= -i^lnj + + + ■'PiJ+l + -ipij-l + h'^i^ij (10) 

■0 = on the boundaries and 

. = -^^il±M ion AD) 

w = -244- (on AB,BC, CD) 
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We will call the polynomial system (fTU]) . pT|) the discrete steady cavity flow 
problem denoted as DSCF{R, s, N). Let the number 2{N — 2)^ of variables in 
DSCF(i?, s, N) corresponding to function evaluations at interior grid points, be 
called the dimension n of the discrete cavity flovi^ problem. Moreover, a 
solution (■0, w)(iV) of the discrete cavity flow problem of discretization N, that 
does not converge to a physical solution of the original continuous cavity flow 
problem for TV — > oo is called a fake solution. 

Several methods have been used to solve the cavity flow problem and the 
steady cavity flow problem numerically (see, e.g., [T], [3], [1], [7] [S], [H], [IH])- 
In this paper, we propose a new method to solve the discrete steady cavity flow 
problem. This method provides solutions sorted by their (discretized) kinetic 
energy 

Remark 1 In we discretize the Jacobian ~ by the central dif- 

ference scheme. It is shown by Arakawa fJ]/ that the central difference scheme is 
the simplest, but the discretized system does not keep important physical invari- 
ants. We study the system DSCF(R, s, N ) as the simplest starting test case. 

Remark 2 We conjecture that the discrete cavity flow problem DSCF(R, s, N ) 
has finite complex solutions. In other words, it defines a zero- dimensional ideal. 
We have checked the conjecture up to N = 5 by Grobner basis computation. 

In the sequel, the boundary velocity s will be denoted by v as long as no con- 
fusion arises with v — —dijj/dx. 

3 Grobner Basis Method 

The Grobner basis method finds all complex solutions of a given system of zero 
dimensional polynomial equations when they are relatively small systems. Be- 
fore discussing the semi-definite programming relaxation method, we will study 
our discrete steady cavity flow problem by the rational univariate representation 
[15], [14], which is a variation of the Grobner basis method. Results will be used 
to tune parameters of the sparse SDP relaxation method. 

The 5x5 mesh is solvable with this method, however the 5x6 mesh case 
is not solvable in one hour by current major implementations (Groebner(Fgb) 
in Maple 11, nd_gr_trace and tolex_gsl in Risa/Asir). The system for the 5x5 
mesh case contains 18 variables and 9 in the 18 appear as linear and the other 
9 as quadratic variables. 

We sort the real solutions by a discretization of the kinetic energy of the 
fluid (mass • velocity^ /2), which is proportional to 




(12) 
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By approximating integral (|12p via the central difference discretization, we ob- 
tain our discrete energy function F, given as 



Fi^P,io) 



E 



2<i,j<N~l 

- E 

2<iJ<N-l 



4/l2 



4/l2 



(13) 



Figure[T]-|3]illustrate the approximations for some velocity vectors (d^p /dy, —dij) / dx) 
of solutions of the discrete steady cavity flow problem obtained by the Grobner 
basis method. We sort real solutions by the discrete energy function. The left 
solution is the minimum energy solution, the center solution is the second en- 
ergy solution, the right solution is the 3rd energy solution. R is the Reynolds 
number, v is the velocity of stream along the boundary, M is the magnification 
factor to display the velocity vector, i.e. we display {Mdip/dy,—Mdip/dx) in 
case M ^ 1. 



Figure 1: i? = 0.01, w = 1, iV = 5, M = l,l/10^1/10^ There are 26 real 
solutions. 




Figure 2: R = 500, v = l,N = 5,M = 2, 1/10, 1/10. There are 14 real solutions. 



It follows from these data that the minimal energy and the energy gap be- 
tween the minimal energy and the second energy seem to decrease when R 
increases. It is interesting that the minimal energy solutions have a vortex of 



5 




Figure 3: i? = 40000, u = 1, = 5, M = 4, 2, 2. There are 20 real solutions. 



clockwise direction, but some of 3rd energy solutions have a vortex of counter- 
clockwise direction, which are apparently fake solutions. 













I 



Figure 4: i? = 500, w = 10, iV = 5, Mag = 2,1/10,1/10. There are 18 real 
solutions. 



4 Sparse SDP Relaxation Method 

The main contribution of this paper is to propose an algorithm that enumerates 
the smallest kinetic energy solutions of the discrete steady cavity flow problem 
DSCF(i?, u, TV) starting with the minimum energy solution. The key element 
of this algorithm is to apply the sparse semidefinite program relaxation method 
(SDPR) to solve the DSCF(i?, v, N). The SDPR for PDEs was proposed in [H] 
and is based on the idea to take the polynomial system derived from a finite 
difference discretization of a differential equation and its boundary conditions 
(for instance: DSCF(i?, w, TV)) as constraints for an optimization problem. After 
choosing a further polynomial function as the objective of the optimization 
problem, a polynomial optimization problem (POP) of the form 

min F{x) 

s.t. gj{x)>0 Vj e {l,...,/c}, (14) 
h,{x)^0 yie{i,...,l}. 
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is obtained. As shown in [T3], polynomial optimization problems derived from 
differential equations satisfy structure sparsity patterns and the sparse SDP 
relaxations due to [5D] can be applied to approximate the solution of POP 
p^ . The crucial point is how to choose the objective function F in POP ^T^ . 
In case that several solutions to a discretized PDE problem exist, the choice 
of the objective function allows to select solutions of particular interest. For 
the cavity flow problem, we are interested in the solution which minimizes the 
kinetic energy p^ . Thus, for the cavity flow problem we yield the previously 
derived function F (fO)) as a canonical choice for the objective function of 

2<i,j<N-l 

(15) 

We define the functions gjj, gf ^ : ^ M as 

Taking into account DSCF(i?, u, N) and the objective function F, we derive the 
polynomial optimization problem, 

min F{ip,Lu) 

s.t. g}^{i;,uj)^0 V2<i,j<iV-l, 
5f,(V',^) = V2<z,j<iV-l, 

^„ =0 y{i,j)e{l,N}x{l,...,N}U{l,...,N}x{l,N}, 
u;i.,=-2^ VjG{l,...,iV}, 
Vj-£{l,...,iV}, 
^..i = -2% Vze {l,...,iV}, 

c.,.^ = -2*^^^-i±^ Vze{l,...,7V}. 

(16) 

We call POP (I16p the steady cavity flow optimization problem CF(R, v, N) 

with Reynold's number R, boundary velocity v and discretization N as param- 
eters. As all polynomials in (|16p are of degree at most two, CF{R,v, N) is a 
quadratic optimization problem (QOP). In fact, a further classification is 
possible for i? = and R 0. 

Proposition 1 a) CF{0,v, N) is a convex quadratic program for any v 

and N . 

b) CF{R, V, N) is non-convex for any v and N, if R^ 0. 
Proof: 

a) In case i? = all constraints are linear. Furthermore, the objective func- 
tion can be written as = ^ ^ ^ F^.^ + Ff ,j , where 

^.(^'-) = (^--^--^) (-2 ^ ) (t!-)- 
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/ 2 -2 \ 

It follows that Fi_j is convex as ( ^ 2 / P*-*^^^^^^ semidefinite with 

eigenvalues and 4. The convexity of Ff^ follows analoguously. Thus, 
F can be written as a sum of convex function and is therefore convex as 
well. The proposition follows. 

b) In case R ^ 0, the equality constraint function g} ^ is indefinite quadratic. 
Thus, CF{R, V, N) is a non-convex quadratic program. □ 

It is our aim to apply the methods proposed in |T3] to solve CF{R, v, N) and the 
underlying discrete steady cavity flow problem, i.e. to approximate the solutions 
of (|16p by solutions of a hierarchy of semidefinite program relaxations SDPR(w) 
constructed in [20' , where w denotes the order of the semidefinite program (SDP) 
relaxation. In theory, the solution of SDPR(w) converges to the optimal solution 
for (fT6|) for w — > oo. Nevertheless, the capacity of present SDP solvers restricts 
the choice of the relaxation order w, as the size of SDPR(i(;) grows rapidly in 
w. However, as pointed out in [20] for many POPs it is suffcient to choose a 
relaxation order w G {wmin, ■ • ■ , + 3} to approximate the POP's minimizer 
accurately. For a general POP p4|) . Wmin denotes the minimal relaxation order, 
which is given by 

J degF rdegg, deg h, 

Wmin = niax < — - — , max — - — , max — - — 

[ 2 l<j<k 2 l<i<l 2 

Therefore, for CF{R,v,N) holds Wmin = 1- 

Remark 3 It is a well known result (c.f. Ill]), that SDPR(l) and {14^ are 
equivalent, in case that the POP {1^-^ is a convex quadratic program. Thus, 
solving CF(0,v, N) is equivalent to solving a SDP. Moreover, it is easy to show 
that the contraints admit only one feasible point when R = 0. 

4.1 Tightening the SDP relaxation and improving the ac- 
curacy 

As stated before, the solution of SDPR(w) converges to the optimizer of the 
POP for w ^ oo. Nevertheless, as the dimension n of CF{R,v, N) is given 
by n = 2(iV — 2)^, choosing a relaxation order w greater than 2 for a medium 
scale discretization N yields a SDP which requires too much memory in order 
to be solved by the used SDP-solver SeDuMi [19]. Therefore, we have to restrict 
ourselves to w = 1, 2 for small scale N, or even to w — I for medium scale N . 
We cannot expect that SDPR(l) or SDPR(2) provide accurate approximations 
to the optimal solution for any R and v. In order to tighten the SDP relaxation 
SDPR(l) and SDPR(2), respectively, we impose lower and upper bounds Ibd''', 
Ibd", ubd*' and ubd'^ G R^' such that 

Ibdf < tpi < ubdf and Ibd^ < lo, < ubd^ V 1 < i < iV^ (^7) 
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holds. 

Furthermore, we may apply additional locally convergent optimization tech- 
niques. For instance Newton's method for nonlinear systems can be applied to 
DSCF(i?, V, N) where the SDPR(it;) solution is taken as the starting point. Or 
alternatively, (|16p is approximated by sequential quadratic programming (SQP) 
[5], again, with the SDPR(u') solution as starting point of the algorithm. Com- 
bining the sparse SDP relaxation with Newton's method or SQP is summarized 
in the scheme: 

Method 1 SDPR method 

1. Choose a boundary velocity v, grid discretization N and Reynolds number 
R'. 

2. Apply SDPR(w) to CF{R' ,v, N) and obtain solution u. 

3. Apply sequential quadratic programming (SQP) to C'F{R' ,v, N) starting 
from u, and obtain u. 

4- Apply Newton's method to DSCF(R' ,v, N ) starting from u or u, and ob- 
tain u. 

We note that Step 3 and Step 4 are optional. 

4.2 Enumeration algorithm for finding the k smallest en- 
ergy solutions 

As mentioned in section 2, we conjecture the number of solutions of the dis- 
crete steady cavity flow problem DSCF(i?, v, A^) is finite, i.e. the feasible set 
of CF{R,v, N) is finite. Method [T] enables us to approximate the global mini- 
mal solution u* = M^^-'* := (^/''^^•'*, w^^-''^) of CF{R,v, N). Beside the minimum 
solution, we are also interested in finding the solution u'-^^* with the second 
smallest kinetic energy, the solution m^'^^* with the third smallest kinetic energy 
or in general the solution uC^)* with the fcth smallest kinetic energy. Based 
on the SDPR method we propose an algorithm that enumerates the k small- 
est kinetic energy solutions of CF{R,v, N). Our algorithm shares the idea of 
separating the feasible set by additional constraints with Branch-and-Bound 
and cutting plane methods that are used for solving mixed integer linear pro- 
grams and general concave optimization problems ^81 . In contrast to the linear 
constraints of those methods we impose quadratic constraints to separate the 
feasible set. Moreover, CF{R,v, N) is a non-convex continuous quadratic op- 
timization problem for R ^ 0. It may be worth investigating in future how 
extensions of Branch-and-Cut methods for certain nonconvex problems [5] can 
be used in our setting. 

Algorithm 1 Enumerate the k smallest solutions: 

Given u'^^~^\ the approximation to the (k — l)th energy solution obtained by 
solving SDPR''-\w). 
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1. Choose ei and 63 > 0. 

2. Choose integers 6^ e {l, . . . , (iV - 2^] . 

3. Add the following quadratic constraints to SDPR'^^^ {w) and denote the 
resulting (tighter) SDP relaxation as SDPR^{w). 

{u, - v^-y := - >4 VI < J < h\, 

K + (Ar-2)^ - "■+(Ar-2)0' K' " ^^'f > 4 VI < J < h^. '^'^> 

4- Solve SDPR^{w) with w = \,2 or larger, if possible. Obtain first approxi- 
mation if^^-f^C^) . 

5. Apply a local optimization technique as for instance Newton's method or 
SQP with it^^^C^) as starting point. Obtain u'^^^ as an approximation to 

6. Iterate steps 1-5. 

The idea of Algorithm [T] is to impose an additional polynomial inequality con- 
straint 

[u, - uY^f := - ^)-^? > el yi<J< bl 

{uj+(N-2)- - "■;(V-2)0' K- - ^j''? > 4 VI < j < bl 

to the POP (fT6|) in iteration k, that excludes the solution u'"'^^ from the feasible 
set of POP (I16p which was obtained in the previous iteration. In case that the 
feasible set of is finite and u''^^ is sufficiently close to u^'^"^)*, the new 
constraint excludes yC^"!)* from the feasible set of and m^'^^* is the new 
global minimizer of (jl6p . Of course, there are various alternatives to step 3 
in Algorithm [l] in order to exclude uC^^^'* from the POP's feasible set. One 
alternative constraint is 

(u, - = {l<i<b), (19) 

where b G {1, . . . , n} , > and Un+i a new additional variable bounded by 
— 1 and 1. It is easy to see that (fT9|) is violated, if u ~ w^^^^)*. Nevertheless, 
it turned out that the numerical performance of (jl9p is inferior to the one of 
for our problem T)SCF{R,v, N). The right tuning of parameters and 6 
is far more difficult for compared to P^ . A second alternative to exclude 
y(fc-i)* g^j.g L-norm constraints such as 




for p > 1. The disadvantage of the constraints (|20p is, they destroy the sparsity 
of the POP p^B]) . as all Ui {i — 1, . . . ,n) occur in the same constraint. Therefore 
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the advantage of the sparse SDP relaxations is lost and the POP can not be 
solved efficiently anymore. An exception is the infinity norm constraint 

II ^ Woo— rnax | — uf^ 1^ (21) 

l<i<n 

In fact, the infinity norm constraints (PT|) are equivalent to the proposed con- 
straints psp for 61 = &2 = -f ■ Because it preserves sparsity and its numerical 
performance is better than the one of we impose (llSp as additional con- 
straints in Algorithm [1] We obtain the following results for Algorithm [TJ 

Proposition 2 Let R, v and N be fixed. Let (w*-^-', . . . be the output 

of the first {k — 1) iterations of Algorithm (Qp. // this output is a sufficiently 
close approximation of the vector of (k — 1) smallest kinetic energy solutions 
(u^^'*, . . . , u'-'^"^)*), and if the feasible set of POP \lb]) is finite, then there exist 
b € {!,..., n} and e G such that the output m'-''') of Algorithm{l\ (for kth 
iteration) satisfies 

u'-''\w) u^^)^ when w ^ 00. (22) 



Proof: As each u'^' is in a neighborhood of u'^'* for all j e {1, . . . ,k — 1}, 
we can choose G {1, . . . , n} and a vector e G M'', s.t. 

Vj e {1, . . . , fc - 1} 3i < 6 s.t. (u, - up'*)' < e,, 

and 

{u, - u[^^*y > e, > fc e {1, . . . , 6} . 

Let CF{R,v,NY''^ denote the CF{R,v,N) with the k systems of additional 
constraints given by step 3 in Algorithm [1] where the fcth constraints are given 
by (HH]) for the constructed b and e. Then it holds 

feas (CF{R, v, Ar)^^)) = feas {CF{R, v, N)) \ ^^u^^^*, u^'^^^)*} . 

Thus, m'^'^)* is the global minimizer of CF{R,v, NY''^ and the global minimum 
is F{u^''^*). As the bounds PT)) guarantee the compactness of the feasible set, 
it holds with the convergence theorem for the sparse SDP relaxations [TU] 

min sdprC^H w) — > min CF{R, V, N)'^''^ = for 

uC^)* for w 00, 

under the assumption that F{u'--^^*) < F{u'^^^*) < . . ., i.e. u'^'^^* is the unique 
minimizer of Ci^(i?,w, AT) □ 

Although we have proven convergence, the capacity of current SDP solvers 
restricts the choice of the relaxation order w to small integers, in our application 
typically to w = 1 or w = 2. Furthermore, we need to choose the parameters e 
and b appropriately, to obtain good approximations of the k smallest kinetic en- 
ergy solutions. In the following numerical examples we will discuss this question 
and show heuristics how to tune these two parameters. 
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4.3 Numerical results 



We will demonstrate the numerical performance of the SDPR(w) and Algorithm 
[1] to enumerate the k smallest solutions. All calculations are conducted on a 
Mac OS X with CPU 2.5GHz and 2 GB Memory. As an implementation of 
the sparse SDP relaxation we use the software SparsePOP [T7] and MATLAB 
optimization toolbox for standard SQP routines in order to improve the accuracy 
of the solution provided by SDPR(i(;). The total accumulated computation time 
in seconds is denoted by tc, the scaled feasibility error of a SDPR solution u' 
w.r.t. the constraints of CF{R, v, N) by Csc- 

4.3.1 CF(4000,1,5) 

In a first setting we choose the discretization 5, i.e. the dimension of 

the POP (|16p is n = 2 ■ 3^ = 18. This dimension is small enough to apply 
the polyhedral homotopy method [H] and the Grobner basis method (Section 
[3]) to determine all complex solutions of DSCF(i?, ti, A^). Therefore, we are 
able to verify whether the solutions provided by Algorithm [T] are optimal. The 
computational results are given in Table [1] Comparing our SDPR results to all 
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2e-10 


4.6e-4 
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1 


le-3 




3 





5 


5e-4 


6.3e-4 




2 


1 


lc-3 




3 





8 


5e-4 


l.Oe-3 





Table 1: Results of Algorithm [T] for Ci^(4000, 1, 5) 



solutions of the polynomial system obtained by polyhedral homotopy method or 
Grobner basis method, it turns out that the solutions u^"^ , u^^' and u^^^ indeed 
coincide with the three smallest energy solutions u'^'^'^* and u'^'^'^* . All three 

solutions are pictured in Figure [5j Note, that the third smallest energy solution 
M^^' shows a vortex in counter-clockwise direction, which may indicate that this 
solution is a fake solution. 

4.3.2 CF(20000,1,7) 

We restrict ourselves to SDPR(l) for solving CF(20000, 1, 7), as the size of the 
SDP relaxation with order w — 2 resulting from this POP of dimension n = 2 ■ 
25^ — 50 is already too large to be solved in reasonable time. The computational 
results for applying Algorithm[T]for different choices of the algorithm parameters 
are enlisted in Table[2] The two parameter settings {e\,h\) = (le — 3, 1) and 
(e},6}) = (le — 6, 5) are not sufficient to obtain an other solution than m^"', 
whereas (ei[,6j) = (le — 5, 5) yields u^^-*, a solution of larger energy. After 
another iteration with {e\,b\) = (le — 5, 5) we obtain a third solution u'^-* of 
even larger energy. The three solutions are pictured in Figure [6l 
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Figure 5: Interior of the solutions u^"-' (left), m^^^ (center) and m^^-' (right) for 
Ci^(4000,l,5) 



k 


w 


4 


^2 


h'l 


"2 


tc 


^SC 


F(wW) 


solution 





1 










2 


3e-7 


3.4e-4 




1 


1 


lc-3 




1 





5 


5e-4 


3.4e-4 




1 


1 


le-6 




5 





5 


6e-6 


3.4e-4 


uW 


1 


1 


le-5 




5 





9 


5e-6 


5.9e-4 




2 


1 


lc-5 




5 





14 


5e-6 


5.2e-3 





Table 2: Results of Algorithm ffl for CF(20000, 1, 7) 



It is interesting to observe that u*^^^ and u*^^^ are one- vortex solutions, whereas 
there seems to be no vortex in the smallest energy solution . 

4.3.3 CF(40000,1,7) 

Next, we examine CF(40000, 1, 7), which is a good example to demonstrate 
that solving DSCF(i?, u, N) and POP is becoming more difficile for larger 
Reynolds numbers. As for the previous problem, the dimension of the POP is 
n = 50, which is too large to be solved by Grobner basis or polyhedral homotopy 
method. Our computational results are reported in Table [3l 



k 


w 


4 


k 

^2 






tc 


^SC 




solution 





1 










3 


2e-7 


3.4e-4 


uW(l) 


1 


1 


5e-6 




5 





7 


6e-9 


7.3e-4 


u(i)(l) 


2 


1 


5e-6 




5 





11 


3c-6 


5.9e-4 


U(2)(l) 


3 


1 


8e-6 




5 





16 


5c-6 


2.3e-4 


U(3)(l) 





2 










5872 


8e-10 


2.6e-4 


u(")(2) 



Table 3: Resuhs of Algorithm [T] for Ci^(40000, 1, 7) 



Solution u'^^^(l) is of smaller energy than ^^-'^^(1), and m("^'(1) is even of 
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Figure 6: Solutions (left), u^^^ (center) and m^^) (right) for CF(20000, 1,7) 



smaller energy than ^^"■'(1). This phenomenon can be explained by the fact, 
that the SDP relaxation with w = 1 is not tight enough to yield a solution 
that converges to u* under the local optimization procedure. The energy of 
m(°)(2) obtained by SDPR(2) is smaller than the one of u(°)(l), but it is not 
the global minimizer as well. In fact. Algorithm [T] with SDPR(l) generates a 
better solution u*^'^^(l) (with smaller energy) in 3 iterations requiring 16 seconds 
computation time, compared to solution u^"^ (2) obtained by applying SDPR(2) 
to CF(40000, 1, 7) requiring 5872 seconds. Thus, applying our enumeration 
algorithm with relaxation order i(j = 1 is far more efficient than the original 
sparse SDP relaxation with w = 2 for approximating the global minimizer of 
POP (fT6|) in this example; our enumeration algorithm fails with the relaxization 
order w = I, but we obtain the global minimizer efficiently by accident. The 5 
different solutions are illustrated in Figure [7] and [51 




Figure 7: Solutions m(°)(1) (left) and ^("^2) (right) 
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Figure 8: Solutions u^^\l) (left), u^^\l) (center) and (right) 



5 Convergence of solutions in case of large Reynolds 
numbers R 

5.1 SDPR(u') for increasing discretization 

In our previous experiments we derived small or even minimal energy solutions 
by Method [1] and Algorithm [1] for various choices of the problem parameters 
R and v with discretization N G {5, 6, 7}. In case that we succeed, applying 
SDPR(w) to CF{R, V, N) yields the minimum kinetic energy solution u* of the 
discrete steady cavity flow problem. The important question arises whether it 
is possible to expand these coarse grid minimum kinetic energy solutions u* to 
finer grids with larger discretization N, i.e. whether these coarse grid solutions 
converge to analytic solutions of the original (continuous) steady cavity flow 
problem for N —>■ oo. As pointed out in, e.g., [4], in case of larger and larger 
Reynolds number R and velocity v the discrete steady cavity flow problem has 
to be solved for finer and finer grids, in order to obtain solutions converging to 
solutions of the continuous problem for N —> oo. In this section we will adress 
the difficult problem to find solutions for CF{R, v, N) converging to continuous 
solutions for large R and pose a related question and a conjecture based on 
computational experiments. As in section 4, the calculations are conducted on 
a Mac OS X with CPU 2.5GHz and 2 GB Memory. 

5.1.1 CF(100, 1,7V) 

We apply SDPR(l) to CF(100, 1, N) and take the solution as starting point for 
Newton's method. Accurate solutions to the discrete steady cavity flow problem 
are obtained for N e {10, 15, 20}. By applying standard grid-refinement meth- 
ods as in [13] , we succeed in extending the solutions to grids of size 30 x 30 and 
40 X 40. The numerical results are enlisted in Table [Hand pictured in Figure 
[9l Thus, it seems reasonable to conclude, that the minimum energy solution 
converges to an analytical solution of the steady cavity flow problem. The dis- 
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Crete steady cavity flow problem has multiple solutions. It is an advantage of 
our method to detect the minimum kinetic energy solution u* converging to an 
analytic solution for N ^ oo. 



N 


w 


^SC 


tc 


F{u*) 


10 


1 


4c- 11 


14 


0.0169 


15 


1 


6e-16 


255 


0.0313 


20 


1 


6e-16 


948 


0.0409 


30 


1 


4e-ll 


1759 


0.0503 


40 


1 


4e-ll 


4156 


0.0554 



Table 4: Results for CF(100, 1, N) for increasing iV 



■ ' / I 



Figure 9: CF(100, 1,7V) solutions for TV = 10 (left), TV = 20 (center), TV = 40 
(right) 



5.1.2 Large Reynolds number R 

In the following we examine CF(10000, 1, TV) for TV e {8, ... , 18}. For all tested 
discretizations we were able to find accurate solutions by SDPR(l) and ad- 
ditionaly applying sequential quadratic programming (SQP). Our results are 
summarized in Table [5] and pictured in Figure [TOl 



TV 


w 




tc 


i^(uW) 


8 


1 


2e-7 


7 


1.5e-6 


10 


1 


3e-10 


21 


3.2e-6 


12 


1 


le-7 


49 


6.0e-6 


14 


1 


5e-9 


99 


l.le-5 


16 


1 


4e-12 


199 


1.9e-5 


18 


1 


2e-8 


501 


3.9c-5 



Table 5: Results for CF(10000, 1, TV) for various TV 
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Figure 10: SDPR(l) solutions of Ci^(10000, 1, N) for iV 8 (left, top), TV = 10 
(center, top), iV = 12 (right, top), TV = 14 (left, bottom), TV = 16 (center, 
bottom) and TV = 18 (right, bottom) 



If we compare the pictures in Figure \10\ it seems the SDPR(l) solution of 
Ci^(10000, 1, TV) evolves into some stream-like solution. Nevertheless, unlike 
the solutions of CF(100, 1, TV), we have not been able to expand this solution to a 
grid of higher resolution by standard interpolation and grid-refinement methods 
so far. It is possible the solution pictured in Figure flOl is a fake solution. 

5.1.3 Large boundary velocity v 

As an example of a setting with larger boundary velocity we study the problem 
CF(500, 10, 7). We apply the SDPR method with relaxation order w — 2 and the 
continuation method, which is a standard method to solve the DSCF we describe 
in section [5T2I and obtain two different solutions for DSCF(500, 10, 7), c.f. Table 
[6] and Figure [TT] It is interesting to observe that both solutions look like stream 
solutions: The SDPR(2) solution with one vortex and the continuation solution 
with two vortices. But as in the previous setting of large R, grid-refinement 
methods to extend these two solutions to higher resolution grids fail. Therefore, 
it seems reasonable to conclude that both solutions are fake solutions. Another 
question is, whether we can derive the continuation solution by Algorithm [TJ As 
for finding the minimum energy solution, we choose w = 2. Choosing 6} > 2 in 
Algorithm [1] generates an SDP that is too large to be solved by the SDP solver 
SeDuMi. The SDP relaxation for 6j = 1 is tractable, but it is too weak to yield 
a solution different than the minimum energy solution for various choices of e\. 

Question 1 Does the minimum kinetic energy solution u* of C'F{R,v, N) con- 
verge to an analytic solution of the steady cavity flow problem for TV — s- 00, even 
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method 


Esc 


tc 




SDPR(2) 


8e-15 


6774 


0.0385 


continuation 


4e-13 




0.0659 



Table 6: Results for Ci^(500, 10, 7) 









\" 1 







Figure 11: SDPR(2) solution (left) and continuation solution (right) for 
CF(500,10,7) 



for large values of R and v ': 



5.1.4 Alternative finite difference scheme 

We mentioned in Remark [1] the simple central finite difference scheme we use 
does not preserve important physical invariants [IJ. Arakawa [IJ proposed an 
alternative finite difference discretization for , that is shown to preserve those 
invariants. We use this alternative scheme to derive an alternative discrete 
steady cavity flow problem ADSCF(i?, w, N) and solve it via the SDPR method. 



In ADSCF(i?, w, A), the finite difference approximation for |^|^ 
is replaced by 



dip du 
dx dy 



m 



- + l ^ ^i+l,j + l) {ipi+lj + V-'ij') 



dy dx Ox dy V-^«' yj > ^ 

+ (Wi+l,j + CJj+lj + l - t^i-lj - UJi^lJ + l) (V'jj + 1 + V'ij) 
■ (Wi+lj--l + t^i+lj — LUi-lJ-l — l^i-l.j) (V'jj' + 1pi,j-l) 



i.j) (V-^j 



-i.i- 



(24) 

It is to be noted that ADSCF(i?, v, N) is less sparse than DSCF(i?, v, N) and 
it is more difficult to derive accurate solutions by SDPR of relaxation order 
1. Nevertheless, we succeed in solving ADSCF(i?, i;, A^) in some instances. For 
example, in Table [7] and Figure [T^ we compare the minimum kinetic energy so- 
lutions obtained for DSCF(5000, 1, N) and ADSCF(5000, 1, A). It is interesting 
that the vortex in the minimum kinetic energy solution for ADSCF(5000,1,A) is 
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preserved for increasing N, whereas the vortex in sohition for DSCF(5000,1,A^) 
seems to deteriorate. 



Problem 


^SC 


tc 


F{u*) 


ADSCF(5000,1,14) 


7e-12 


1304 


1.8e-4 


ADSCF(5000,1,16) 


5e-10 


2802 


3.1e-4 


DSCF(5000,1,14) 


le-11 


419 


5.6e-4 


DSCF(5000,1,16) 


3c-10 


768 


l.le-4 



Table 7: Results for solving ADSCF(5000,1,N) compared to DSCF(5000,1,N) 




Figure 12: Solutions for ADSCF(5000,1,14) (top left), ADSCF(5000,1,16) (top 
right), DCSF(5000,1,14) (bottom left) and DCSF(5000,1,16) (bottom right) 



5.2 Solutions of CF{R,v, N) for increasing Reynolds num- 
ber R and velocity v 

For small Reynolds numbers, we have seen that the minimum kinetic energy so- 
lution converges to an analytic solution for ^ cx) by applying grid-refinement 
methods. In order to adress Question [1] and to understand why convergence 
to the analytic solution is a lot more difficult to obtain for large R and v, we 
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examine the behavior of the minimum energy solution of the polynomial system 
DSCF{R,v, N) and CF{R,v, N), respectively, for increasing Reynolds number 
R. 

5.2.1 Minimum kinetic energy solution for increasing R 

To solve the discrete steady cavity flow problem, the proposed SDPR method, 
Method [H is one possibility to find an appropriate starting point for Newton's 
method. If w is chosen sufficiently large, the output u of Method[l]is guaranteed 
to accurately approximate the minimum energy solution u* of CF(R' , v, N) and 
DSCF(i?', u, iV), respectively. In order to show the advantage of the SDPR 
method we compare our results to solutions of DSCF(i?', v, N) obtained by 
a second procedure: In case R = the discrete steady cavity flow problem 
DSCF(_R, V, N) is a system of linear equalities, which has an unique solution 
uo{v,N), or short uq. Beside solving the linear system, one way to obtain this 
solution is solving CF(0, v, N), which is equivalent to solving an SDP as pointed 
out in Remark O 

Method 2 Naive homotopy-like continuation method 

1. Choose a boundary velocity v, grid discretization N, Reynolds number R' 
and step size AR. 

2. Solve DSCF(0,v, N ), i.e. a linear system, and obtain u". 

3. Increase R^-^ by AR: R'' ^ R^-^ + AR 

4. Apply Newton's method to DSCF(R^ ,v,N ) starting from u''^^ . Obtain 
solution u'' as an approximation to the discrete cavity flow problem's so- 
lution. 

5. Rerate 3. and 4- until the desired Reynold's number R' is reached. 

We call MethodlHthe continuation method. In fact, it is one of the standard 
methods to find a solution for the steady cavity flow problem. Note, the contin- 
uation method does not necessarily yield the minimum kinetic energy solution of 
DSCF(i?, V, N). In all numerical experiments the boundary velocity v is fixed to 
V — 1. Let u*{R, N) denote the global minimizer of CF{R, 1, N), the minimum 
energy £',„i„(i?, iV) is given by EminiR,N) F{u*{R,N)). Obviously, it holds 
-E,„i„(0, N) = F{uoiN)), Figurelllshows £;,„i„(0, N) for N ranging from 5 to 20. 

In a next step the solution of DSCF(i?, v, N) obtained by the continuation 
method starting from uq is denoted as u{R), and its energy as Ec{R,N) := 
F{u{R,N)). As illustrated for = 5 and = 7 in Figure [Til it is possible for 
all R to find a continuation t2 of mq. For N = 5 the dimension n of the discrete 
steady cavity flow problem is n = 18. This dimension is small enough to solve 
a polynomial system by Grobner basis or polyhedral homotopy method and to 
determine all complex solutions of the system. Therefore, we can verify whether 
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Figure 13: E^in{0,N) for N e [5,20] 



SDPR combined with SQP detects the global minimizer of CF{R, v, N) or not. 
It is worth pointing out, that we are able to find the minimum energy solution of 
the CF{R, V, N) by applying the SDP relaxation method, whereas this solution 
cannot be obtained by the standard continuation method. We observe SDPR(l) 
is sufficient to detect the global optimizer for R < 10000, and for R > 20000 
the global optimizer is obtained by SDPR(2), which is reported in Tabled 




Figure 14: Ec{R) and E^in{R) in case N = 5 (left) and = 6 (right) 



In case of iV = 6 and N = 7 the dimension of the polynomial system is 
too large to be solved by Grobner basis or polyhedral homotopy method for 
i? > 0. For iV = 6 the continuation method, SDPR(l) and SDPR(2) yield the 
same solution for all tested R as pictured in Figure fT4l (right). And in case of 
N = 7 the continuation solution u{R) is detected by SDPR(l) as well, except 
the case R = 6000, where a solution with slightly smaller energy is detected by 
SDPR(l), as documented in Tableland illustrated in Figure [TSl 

Summarizing these results, F{uo{R, N)) > F{u{R, N)) for any of the tested 
i? > 0. It is an advantage of our approach to show, u{R, N) is in general 
not the optimizer of CF{R, 1, TV) for increasing R. In fact, for some settings we 
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R 


Nc 


Nr 


Ec{R) 


EsDPR(l) 


EsDPR(2) 







1 


1 


0.0096 


0.0096 


0.0096 


0.0096 


1 


35 


23 


0.0096 


0.0096 


0.0096 


0.0096 


10 


37 


17 


0.0094 


0.0094 


0.0094 


0.0094 


100 


37 


13 


0.0030 


0.0030 


0.0030 


0.0030 


200 


37 


11 


0.0013 


0.0013 


0.0013 


0.0013 


500 


37 


13 


6.2e-4 


6.2e-4 


6.2e-4 


6.2e-4 


1000 


37 


13 


5.4e-4 


5e-4 


5e-4 


5e-4 


2000 


37 


13 


6.2e-4 


6.2e-4 


6.2e-4 


6.2e-4 


3000 


38 


18 


6.5e-4 


4.8e-4 


4.8e-4 


4.8e-4 


4000 


37 


17 


6.3e-4 


4.6e-4 


4.6e-4 


4.6e-4 


6000 


36 


16 


5.7e-4 


4.5e-4 


4.5e-4 


4.5e-4 


8000 


36 


16 


5.2e-4 


4.5e-4 


4.5e-4 


4.5e-4 


10000 


35 


17 


4.7e-4 


4.5e-4 


4.5e-4 


4.5e-4 


20000 


35 


17 


4.5e-4 


4.5e-4 


3.3e-4 


3.3e-4 


30000 


35 


17 


4.5e-4 


4.5e-4 


2.5e-4 


2.5e-4 


50000 


35 


17 


4.5e-4 


4.5e-4 


1.7e-4 


1.7e-4 


70000 


35 


16 


4.5e-4 


4.5e-4 


1.2e-4 


1.2e-4 


100000 


34 


16 


4.5C-4 


4.5e-4 


8.8e-5 


8.8e-5 



Table 8: Numerical results for CF{R, 1, 5) 



R 





50 


100 


500 


2000 


4000 


6000 


8000 


10000 


Ec{R) 


2.0e-2 


1.4e-2 


7.7e-3 


9.3e-4 


4.5e-4 


4.1e-4 


3.7e-4 


3.5e-4 


3.4e-4 


-E'SDPR(l) 


2.0e-2 


1.4e-2 


7.7e-3 


9.3e-4 


4.5e-4 


4.1e-4 


3.6e-4 


3.5e-4 


3.4e-4 



Table 9: Numerical results for CF{R, 1, 7) 



obtain far better approximations to the minimum energy solution than N). 
Furthermore, Emin{R) and Ec{R) are both decreasing in R. The behavior of 
Ec^ EsDPR and -Emin coincides for all chosen discretization N and motivates 
the following conjecture. 

Conjecture 1 Let boundary velocity v and discretization N be fixed. 

a) F{uo{v,N))=E^in{0,v,N)>E„,in{R,v,N)>0 VR>0. 

b) E^in{R,v,N) for R^ oo. 

As an application, Conjecture [1] can be used as a certificate for the non- 
optimality of a feasible solution u' of CF{R,v, N) in case F{u' {R,v, N)) > 
-Emin(0, w, iV). In fact, as it seems to be always possible to extend uq via con- 
tinuation method, u{R,v,N) can serve as a non-optimality certificate in case 
F{u'{R,v,N)) > F{u{R,v,N)). 
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Figure 15: Ec{R) and Esdpr{i){R) in case N — 7 



5.2.2 Stability analysis 

Finally, we examine the stability of the minimum kinetic energy solution. We fix 
N and v — 1 and increase the Reynolds number R. Let J{u) denote the Jacobi 
matrix of the polynomial system T)SCF{R,v, N) at the solution w, A,nax(u) its 
maximal eigenvalue and its number of positive eigenvalues. A solution u 
to DSCF{R,v, N) is called stable, if all eigenvalues of J{u) are non-negative, 
otherwise it is called unstable. In case = 10 and N = 20, we observed that 
the minimum kinetic energy solution u*{R) is stable for small _R, and as R is 
increased and exceeds some threshold R'{v,N), u*{R) becomes unstable. See 
Table [TO] and Figure [TOl 



R 


N 


^scaled 


tc 


Amax(w) 




100 


20 


9e-16 


491 


-0.0636 





750 


20 


5e-14 


366 


-0.0615 





775 


20 


6e-16 


486 


-0.0014 





776 


20 


4e-16 


486 


0.0010 


2 


800 


20 


7e-16 


486 


0.0580 


2 


1000 


20 


lc-12 


527 


0.5124 


2 


300 


10 


6e-16 


10 


-0.3069 





350 


10 


6e-16 


11 


-0.1360 





400 


10 


3e-16 


9 


0.1217 


2 


500 


10 


4c-16 


9 


0.6043 


4 



Table 10: The stable solution u*{R) becomes unstable for increasing R. 
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Figure 16: u*{R) for R = 100 (left), R ^ 776 (center) and R = 800 (right) 



6 Conclusion 

We proposed an algorithm to enumerate all solutions of the discrete cavity 
flow problem with respect to their kinetic energy. Our algorithm takes advan- 
tage of the sparse semidefinite relaxation method (SDPR) in order to find a 
good starting point for Newton's method or sequential quadratic programming. 
We can guarantee the convergence of the algorithm's output to the smallest 
kinetic energy solutions of the polynomial system, if the order of the SDPR 
tends to infinity. Our numerical experiments for various choices of R and v 
have demonstrated that it is sufficient to apply SDPR of order one or two, to 
succeed in obtaining accurate approximations to the smallest energy solutions 
of the discrete cavity flow problem by our enumeration algorithm. In case of 
small Reynolds numbers our algorithm allowed another interesting observation: 
Among all solutions of the polynomial system given by the discrete cavity flow 
problem, the minimal kinetic energy solution converges to an analytic solution 
of the continuous steady cavity flow problem. In case of large Reynolds num- 
ber R we are not able to extend our coarse grid solutions to a finer grid, yet, 
although many of them look like stream solutions when the kinetic energy is 
small. It is known that the set of solutions of the discrete cavity flow problem 
contains lots of non-physical solutions or fake solutions, but there has been no 
systematic study of the discrete cavity flow problem as a polynomial system so 
far. Moreover, the more interesting stream-like solutions of the discrete steady 
cavity flow problem are usually among the 3rd or 4th smallest kinetic energy 
solutions. Our enumeration algorithm based on the SDPR method provides a 
powerful tool to detect the smallest energy solutions one by one, which is a 
strong advantage compared to the existing methods. The further analysis of 
the polynomial system derived from the steady cavity flow problem for large 
Reynolds number R will remain an interesting topic in future. Also, the conjec- 
ture that the mininum kinetic energy converges to zero for increasing R is left 
for future research and may constitute an interesting property of the minimum 
energy solution, which does not converge to a zero solution itself. 

To conclude, we think that the polynomial system of the discrete steady 
cavity flow problem is challenging for the community of solvers of polynomial 
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systems and numerical algebra. Another interesting challenge is to solve the 
discrete steady cavity flow problem derived by the alternative finite difTerence 
discretization of the Jacobian proposed by Arakawa [T] . Our first computational 
results suggest it is worth further pursuing this alternative. For its observed and 
described properties the discrete steady cavity flow problem will be a good test 
problem to validate new techniques for solving systems of algebraic equations 
and inequalities. Furthermore, as solving the cavity flow problem for large 
Reynolds numbers R and velocities v remains an active field of research, we 
believe that our numerical results may be instructive for audiences in the com- 
munity of numerical analysis for fluid dynamics to understand fake solutions in 
partial differential equations. 
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